/*
 * GaseousActivity.cpp
 *
 *  Created on: 20 Aug 2011
 *      Author: Allan
 */

#include "GaseousActivity.h"

// GeoReact includes
#include "GaseousModels/SpycherReedCO2Model.h"

GaseousActivity::GaseousActivity() 
{}

GaseousActivity::GaseousActivity(const vector<string>& species) :
AuxiliaryActivity(species, IdealGaseousModel())
{
	// Set the Spycher and Reed (1988) model for CO2(g) (only effective if CO2(g) was provided)
	SetActivityCoefficient("CO2(g)", SpycherReedCO2Model());
	
	// Set the Spycher and Reed (1988) model for H2O(g) (only effective if H2O(g) was provided)
	SetActivityCoefficient("H2O(g)", SpycherReedH2OModel());
}

GaseousActivity::~GaseousActivity()
{}

const VectorXd GaseousActivity::Activities(double T, double P, const VectorXd& n) const
{
	VectorXd a(n.rows());
	
	// The reciprocal of the total number of moles of the gaseous species
	const double ntotalRec = 1.0/n.sum();
	
	// Calculate the activity of each gaseous species as the product of its molar fraction with its activity coefficient and pressure
	for(int i = 0; i < a.rows(); ++i) 
	{
		// The molar fraction of the i-th gas
		const double xi = n[i] * ntotalRec;
		
		a[i] = P * xi * activitycoeffs[i](T, P, n);
	}
	
	return a;
}

const ActivityCoefficient IdealGaseousModel()
{
	// The ideal model for the gaseous activity coefficient (gamma = 1.0)
	auto ideal = [](double T, double P, const VectorXd& n) { return 1.0; };
	
	return ideal;
}

const ActivityCoefficient SpycherReedCO2Model()
{
	return ActivityCoefficient(SpycherReedFugacityCoefficientCO2);
}

const ActivityCoefficient SpycherReedH2OModel()
{
	return ActivityCoefficient(SpycherReedFugacityCoefficientH2O);
}
